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We use a lattice Boltzmann algorithm for liquid-gas coexistence to investigate the 
steady state interface profile of a droplet held between two shearing walls. The 
algorithm solves the hydrodynamic equations of motion for the system. Partial 
wetting at the walls is implemented to agree with Cahn theory. This allows us to 
investigate the processes which lead to the motion of the three-phase contact line. 
We confirm that the profiles are a function of the capillary number and a finite size 
analysis shows the emergence of a dynamic contact angle, which can be defined in 
a region where the interfacial curvature tends to zero. 
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1. Introduction 

When a liquid drop in equilibrium with its vapour is placed in contact with a 
flat surface the equilibrium shape defines a static contact angle, 0^,, between the 
liquid-vapour interface and the surface. 6^ is determined by a balance between the 
fluid-fluid and solid-fluid interactions. The surface may be wet [6^ = 0), partially 
wet (0 < < tt) or dry (6'^, = tt) (de Gennes 1985). 

If the droplet is pushed over the substrate its steady state proflle defines two 
dynamic contact angles (advancing and receding) different, in general, from the 
static angle. The shape of the moving droplet is difficult to investigate analytically 
because the classical continuum hydrodynamic equations of motion with the usual 
no slip condition at the surface predict a singularity in the stress at the contact line 
(Huh & Sriven 1971; Dussan & Davis 1974). 

In analytic work to remedy this problem it has been proposed that a slip con- 
dition may hold at the three phase point and matched asymptotic expansions have 
been used to describe the interfacial configuration and flow fields near to the bound- 
ary (Hocking 1977). More recently Seppecher (1996) used a diffuse interface model 
with a no slip condition to provide a mesoscopic description of contact line motion. 
By matching a diffuse interface solution close to the contact line to a sharp interface 
solution far from it, Seppecher was able to show how the dynamic contact angle is 
related to the capillary number and to the static angle. 

Numerical modelling of contact line motion has also proved difficult because 
of the widely differing length scales involved in the problem. Molecular dynamics 
simulations give useful information on the local boundary conditions but are unable 
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to reach length and time scales on which the dynamic contact angle can be mea- 
sured. Numerical solutions of the Navier-Stokes equations which asssume a sharp 
interface again suffer from the problem of an infinite stress at the contact line. 

It therefore seems appropriate to investigate the extent to which emergent 
mesoscale modelling techniques allow modelling of the dynamics of droplet mo- 
tion. These appproaches incorporate a diffuse interface, thus removing any contact 
line singularity. They are also able to address longer length and time scales than 
molecular dynamics, albeit at the expense of a lack of molecular detail. A main aim 
of this paper is to ask whether it is possible to meaningfully define an advancing and 
receding dynamic contact angle within the framework of one mesoscale approach, 
lattice Boltzmann simulation (Chen & Doolen 1998; Sued 2001). 

We use the one-component, non- ideal lattice Boltzmann model proposed by 
Swift et al (1996) to investigate motion of the contact line. The simulations model 
a droplet held between parallel plates, which impose a shear flow on the system. 
We have developed boundary conditions which allow us to control the wetting 
behaviour at the plates. In particular, the static contact angle agrees with that 
predicted analytically by Cahn theory (Cahn 1977). The simulations show that, 
in the steady state, a dynamic contact angle can be defined away from the wall. 
However, even in two dimensions, rather large lattices are needed to do this. 

This paper is organised as follows. In §H we outline the lattice Boltzmann scheme 
used, ^ describes the boundary conditions necessary to simulate wetting and §|| 
presents the results of the simulations. Finally in §^ we discuss the results and draw 
conclusions. 



Lattice Boltzmann simulations (Chen & Doolen 1998; Succi 2001; Swift et al 1995, 
1996) solve the Navier-Stokes equations by following the evolution of a set of distri- 
bution functions, /o-i(x, t), which represent the density at time t and lattice site x 
which is travelling with velocity e^i. The velocity vectors eo-i are such that the dis- 
tribution functions advect to neighbouring lattice sites x -I- Ax in the time interval 
Af. The velocity vector subscripts a, i and a are used to specify a vector's mag- 
nitude, direction and Cartesian components. For a square lattice and a nine-speed 
lattice Boltzmann model which we consider here, a — Q for the zero velocity vector, 
CT = 1, i = 1, 2 . . . 4 for vectors to nearest- neighbour sites and f7 = 2, i = l,2...4 for 
vectors to next-nearest neighbour sites. Physical quantities are related to moments 
of fai- The density n and fluid velocity u are defined by 



The distributions fcri are evolved according to a lattice Boltzmann equation 
assuming a single relaxation time approximation 



2. Method 




(2.1) 




(2.2) 



/<,,(x + e<,,At, t + M)^ /^,(x, t) 




(2.3) 
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where r is the relaxation time and f^j is an equiUbrium distribution function. The 
equiUbrium distribution determines the physics inherent in the simulation. A power 
series in the local velocity is assumed (Swift et al 1996) 

where summation over repeated Cartesian indices is understood. 

The coefficients A^, B„, C^, D„ and Gaap are determined by placing constraints 



on the moments of /^f . In order that the collision term in equation (2.3) conserves 
mass and momentum the first two moments of are constrained by 

E/-="' (2-5) 

ai 

rji<icia = nua- (2.6) 

ai 

The next moment of is chosen such that the continuum macroscopic equa- 



tions approximated by the evolution scheme (2.3) correctly describe the hydrody- 



namics of a one-component, non-ideal fluid. This gives 

fl1<iaiaec,if} = Pap + nUaUp + v[uadp{n) + upda{n) + u^d^{n)5ap\ (2.7) 

ai 

where v = {t — At/2)/3 is the kinematic shear viscosity and Pap is the pres sure 
tensor. The first formulation of the model omitted the third term in equation ( ^.7[) 
and was not Galilean invariant. Holdych et al (1998) showed that the addition of 
this term led to any non-Galilean invariant terms being of the same order as finite 
lattice corrections to the Navier-Stokes equations. In order to fully constrain the 
coefficients A^r, B^, C^, and Gaap a fourth condition is needed (Huo et al 1995) 

fale^iaCaipeai^ = ^{Ua5fj-f ^ UfjSaj + U^5af}). (2.8) 



The analysis of Holdych et al (1998) shows that the evolution scheme (O) approx- 
imates the following Navier-Stokes level equation: 

dtijlUa) + dpinUaUp) = -dpPap + vdf3[n{di3{Ua) + da{up) + d~f{u^)Saf3)] 
— ivdp \uad^Pp^ + UpdryPa-y + dj{nUaUf3Uj)'\ 

- ivdp [{dnPap)dy{nUy)\ 
- iv"^ dp[uady{updy{n) + Uydfs{n) -f uxdx{n)5ap)\ 
— 'iu'^di3[ui3dy{uady{n) + u^da{n) + u\d\{n)5ai3)\ 

+ 3^/^9/3 [dt{uadi3{n) + upda{n) + uxdx{n)5ap)] ■ (2.9) 

The thermodynamics of the fluid enter the lattice Boltzmann simulation via the 
pressure tensor P^p (Rowlinson & Widom 1982). For a system without surfaces 
the equilibrium properties of the fluid can be described by a Landau free energy 
functional of the form (Landau & Lifshitz 1958) 

= j dV[ij{T,n) + ^{danf] (2.10) 
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where k is related to the surface tension, and (Rowhnson & Widom 1982) 

^P{T, n) = W{T, n) + ^ii,{T)n - pb{T). (2.11) 

Here, fib and pi, are the chemical potential and pressure in the bulk, and W is a 
non-negative function of n that vanishes, along with dW/dn, when n is equal to the 
liquid bulk density ti/ or to the gas bulk density Ug. It then follows that (Rowlinson 
& Widom 1982) 

Pa/3(x) = Sapp{x) + K.{dan){df}n) (2.12) 

with 

p[x.) = pq — Knd^-f-n — '^{d^n)'^ (2.13) 

where po = ndn^'iT, n) — ip{T, n) is the equation of state of the fluid. However, to 
incorporate wetting the choice of free energy must be generalised to include surface 
terms and we now consider this case. 



3. Wetting 

When a liquid-gas interface meets a solid wall the angle, 9^, between the interface 
and the wall, measured in the liquid, is determined by the liquid-gas, solid-liquid 
and solid-gas surface tensions, tr, agi and agg according to Young's equation (Young 
1805) 

cosOw = — . (3.1) 

(T 

In this section our aim is to define lattice Boltzmann boundary conditions which 
reproduce Young's equation in equilibrium. The solid-gas and solid-liquid surface 
tensions will be related to an additional term in the Landau free energy functional 
which describes the interactions at the surface between the solid and the fluid. To 
this end we follow Cahn (1977) and introduce an additional surface term into the 
free energy 



^s = JdS $(n,) (3.2) 



where is the fluid density at the wall. Following de Gennes (1985) we expand 
^{us) as a power series in Ug. In addition, and following Seppecher (1996), we keep 
only the first order term and write $ = —(piUg since this turns out to be sufficient 
for the partial wetting scenarios that we will want to consider. 

We now summarise how the wetting angle 9^ can be written in terms of the 
wetting potential 4)i (Cahn 1977; Papatzacos 2001). To do this we calculate the 
surface tensions by considering a one dimensional problem where one phase of a 
non-ideal fluid occupies the region x > with a solid wall at a; = 0. Far from the 
wall the fluid density will be rif, (where rii, is the bulk liquid density ni or the bulk 
gas density rig), while at the wall n = U s with Ug unde term ined as yet. In one 



dimension the free energy functional ( p.lC| ) together with (3^) reduces to 



+ = / dx[V'(T,n)-|-|(9,n)2] (3.3) 
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Minimising (3.3) subject to natural boundary conditions leads to two conditions 



dip (Pn 

— - /i6 - 1^-^ ==0 for a; > 0, 



/dn^ d$(ns 
^dx dn. 



at a; = 0. 



(3.4) 
(3.5) 



Equation ( p.4D is the usual Euler-Lagrange equation and (^) is a boundary con- 
dition valid at .T = 0. A first integral for equation (3.4) is 



2^dx' l^bn + Pb^Win) 



(3.6) 



suppressing for now the T dependence of W . We can determine by substituting 
(3.6) into (^) giving 



-(t)i = ±^/2KW{ns). 



(3.7) 



Consider first (pi > 0- Equation (3.7) has four solutions (rii < n2 < < 114) if 
(pi is smaller than the height of the double well function defined by ^/2kW (Cahn 
1977). The value of Us is obtained from one of these four solutions as the one which 
minimises the solid-fluid surface tension 



-(pins 



V2kW dr 



(3.8) 



For (pi small enough it can be shown that the minimising solutions are n2 if nf, = Ug 
and 714 if = ni. This gives expressions for the suface tensions 



-<pin2 ■ 

-01714 ■ 



/ V2kW dn, 

J Ug 

/ V2kW dn. 

J ni 



(3.9) 
(3.10) 



Similarly, ii <pi < the minimising solutions are rti if rif, = Ug and 77.3 if 7ib = ni 
and the solid-fluid surface tensions are 



'sg 



-(pins 



m 



V2kW dn, 
V2kW dn. 



(3.11) 
(3.12) 

(3.13) 



To find a closed form for 9^,; we now choose the excess free energy function 
W{n,T) to be 



The liquid-gas surface tension a follows from (^^) as 

V2kW dn. 



(3.14) 
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where v = (n — n^jric and r = (Tj, — T)/Tc. Tc, Pc and ric are the critical tem- 
perature, pressure and density respectively, and /3 is a constant. With this form for 
W{n,T) Young's law reduces to 



where H, = Pt^2k,Pc. 

In order to implement this scheme in a lattice Boltzmann simulation equation 
( ^.5| ) is imposed on lattice sites which represent the wall. Our approach in similar to 
that recently taken by Desplat et al ( 2001) in introducing wetting boundary condi- 



tions for a binary mixture. Since (3^) is an equilibrium condition, it is appropriate 
to impose it through the equilibrium distribution function, /^*. The coefficients of 
/^f depend on the local values of n, Vn and V^n which, in the bulk, are calculated 
using standard finite difference methods. For a wall parallel to a lattice direction, 
it is the perpendicular components of Vn and V^n which must be calculated using 



equation (3.5). For dxn (where x is the perpendicular direction to the wall) we 
use the calculated value of For (9^n, we use the standard right-handed finite 

difference formula 

dlnl.. - (3.16) 

where n[ is dxn at lattice site i. In this formula we substitute for tiq using equation 
( ^.5| ) and calculate n'l using a standard centred finite difference formula. Finally we 
have found empirically that the best choice for n'2 is a left-handed finite difference 
formula taken hack into the wall 

n^K ^ ■ (3.17) 

Using this scheme to evaluate /^f at wall sites it is possible to control the wetting 
angle at any flat wall. To validate the method we have simulated a droplet of liquid 
in equilibrium with its gas on a solid surface. Figure |l] shows how the observed 
contact angle varies with O. The agreement between simulation and theory is good 
and therefore we take this method as a basis for simulations of spreading. 



4. Contact line dynamics 

We now present the results of simulations exploring the motion of a contact line. 
Our aim is to show how a dynamic contact angle can be measured using the lattice 
Boltzmann method. 

In his analytic work on contact line motion, Seppecher (1996) considered a 
droplet on a surface and defined three regions of flow for his matched solution. 
In the inner zone, a semicircle of radius Seppecher used the Cahn-Hilliard 
equation to solve for the flow field with the condition that the interface intersected 
the wall at the static contact angle, S^. Interfacial curvature is concentrated in 
this region and is the cause of mass transfer across the interface, which enables 
the interface (but not the fluid) to slip relative to the surface. Seppecher matched 
this inner solution to a solution of the Stokes equation in the intermediate zone. 
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In the intermediate zone, an annulus of outer radius -R2(3> the interface was 
considered to be flat and have zero width. The angle this part of the interface makes 
with the surface is taken to be the apparent dynamic contact angle, $ — the angle 
observable in experiments. The intermediate zone is matched to the outer zone 
which is the region far from the contact line. In the external zone the curvature 
radius of the interface is much larger than i?2 and the external forces on the droplet 
and geometry of the whole domain determine the flow. 

Using diffuse interface simulations, and a droplet of a computationally feasible 
size, the intermediate zone, which defines the dynamic contact angle, is masked 
by the curvature of the droplet. Therefore in our simulations we have employed a 
different geometry which allows us to probe the intermediate zone within the limits 
of the computational resources available. 

Our simulations model a two dimensional droplet held between two partially 
wetting walls at a; = and x ~ L^ — l, where is the number of lattice points in the 
X direction. The wetting properties of the two walls may be varied independently. 
The system is Ly lattice sites in the y direction where periodic boundary conditions 
close the system. The droplet is allowed to reach equilibrium with its gas, and then 
a shear flow is applied, by imposing velocities ±Vo on the walls. The system achieves 
a steady state and the interface profile is recorded by measuring the angle, 0{x), 
the interface makes with the wall at x = (see figure ||). 

Initial simulations focussed on confirming the dependence of the interface profile 
on the capillary number, Ca = rjU/a where rj is the viscosity, U is the speed of 
the interface relative to the wall and a is the surface tension. For small systems 
{Lx = Ly = 100) we varied the capillary number from a reference value by inde- 
pendently varying the viscosity (case A) or the shear rate (case B): The aim was 
to show that 9{x) depends only on the product i]Vo. Figure || shows steady state 
interface profiles for a reference value of the capillary number, Ca = 0.085, and for 
capillary numbers four and seven times greater than this value. This figure reveals 
that 0{x) scales approximately with capillary number, although there are some sur- 
prising discrepancies. However, the discrepancies between cases A and B for each 
capillary number can be explained as being due to spurious velocities in the lattice 
Boltzmann model (which are a lattice artifact). The case of increasing r (case A) 
converges on the results of r = 0.8 (case B) if the resolution of the grid is increased 
(Briant, unpublished data). This value of r minimises the spurious velocities due 
to cancellation of higher order terms in the Chapman-Enskog expansion, as noted 
by Swift et al (1996). 

We now address how to measure a dynamic angle. From Seppecher (1996), 
we expect an intermediate zone, where the curvature of the interface is zero, to 
develop far from the contact line. Therefore we perform the shearing experiments 
described above and vary the number of lattice sites in the x direction, to see 
whether such an intermediate zone can be observed. Figure ^ shows the interface 
profiles achieved for two neutrally wetting walls {6^ = 90°) and lattice sizes of 
Lx = 150, 200, 250, 275, 300, 350 and 400. 

Figure ^ shows that as the distance between the walls is increased, the angle in 
the central region of the domain approaches a limiting value. This is more clearly 
indicated in figure where we have plotted 6{x/ L^ = 0.5) against L~^ for each 
set of data. Figure |g shows that the midpoint angle is converging to a well-defined 
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value for an infinite system. Identifying the extrapolated value as the dynamic angle 
gives $ = 110°. 

5. Discussion 

In this paper wc have described a lattice Boltzmann scheme for the simulation of 
partial wetting and of moving contact lines in a liquid gas system. Wetting boundary 
conditions have been defined which enable the contact angle of the interface to be 
controlled in a way consistent with Cahn theory. Wc have used results for the shape 
of a droplet held between sheared parallel plates, together with a finite size analysis, 
to investigate the dynamic contact angle. The approach demonstrates that a lattice 
Boltzmann model with hydrodynamic no slip boimdary conditions can lead to a 
well defined dynamic contact angle. However, the limitation of a diffuse interface 
approach is that it is difficult to model large domains because of the magnitude 
of the computational demand. Therefore measuring a meaningful dynamic contact 
angle for, say, a droplet moving along a surface would prove very difficult. 

Future directions for the work described in this paper include the study of dif- 
ferent wetting conditions at the sheared walls and finding the functional relation- 
ship between the capillary number and the dynamic contact angle. It is of interest 
that recent molecular dynamics simulations of binary fluids (Denniston & Rob- 
bins 2001) have questioned the validity of no slip boundary conditions. Although 
mesoscale simulations are unable to investigate the molecular basis for imposing 
a given boundary condition they provide an ideal tool for understanding how a 
change in the boundary conditions employed will affect the interface profile. 

Jones et al (1999) have used a different mesoscale approach, dissipative particle 
dynamics, to investigate how a droplet in a binary fluid is pulled from a wall under 
shear flow. Dissipative particle dynamics simulations differ from lattice Boltzmann 
in that they include fluctuations. It would be interesting to compare results from the 
two approaches for the geometry described here. Recently Fan et al (2001) have used 
the two component lattice Boltzmann method developed by Shan & Chen (1993) 
to study the contact line dynamics of a pressure driven binary fluid in a capillary 
tube. They found that the cosine of the contact angle varies linearly with the speed 
of the contact point for two different wetting conditions. Other lattice Boltzmann 
investigations of contact line motion for binary fluids have been given by Grubert 
& Yeomans (1999) and Desplat et al (2001) using a free energy approach similar 
to that given here. These authors point out that for binary systems interspecies 
diffusion is important in facilitating the motion of the interface near the wall. 

We thank A. Balazs, C. Denniston, I. Pagonabarraga and P. Warren for helpful discussions. 
AB acknowledges EPSRC research studentship 99315535 and a CASE award from Unilever 
pic. 
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Figure 1. Static wetting angle 6^ plotted as a fu nctio n of dimensionless wetting potential 
n. Curve: theoretical relation, equation ( ^.15| ); circles; simulation results. 
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Figure 2. Sheared interface profiles showing the definition of 9{xo)- 
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Figure 3. Interface profiles d{x) plotted against :c/Lx for a reference capillary number, 
Ca = 0.085 (A), and capillary numbers four (case A:*; case 3:0) and seven (case A:o; 
case B:o) times greater. 
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Figure 5. Mid point interface angles 9{x/Lx = 0.5) plotted against l/L^- 
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